Negative heat capacity of sodium clusters 



Juan A. Reyes-Nava, Ignacio L. Garzon, and Karo Michaelian 
Instituto de Fisica, Universidad Nacional Autonoma de Mexico, 
Apartado Postal 20-364, 01000 Mexico D.F., Mexico 
(Dated: February 2, 2008) 

Heat capacities of Najv, N = 13, 20, 55, 135, 142, and 147, clusters have been investigated using 
a many-body Gupta potential and microcanonical molecular dynamics simulations. Negative heat 
capacities around the cluster melting-like transition have been obtained for N = 135, 142, and 147, 
but the smaller clusters (N = 13, 20, and 55) do not show this peculiarity. By performing a survey of 
the cluster potential energy landscape (PEL), it is found that the width of the distribution function 
of the kinetic energy and the spread of the distribution of potential energy minima (isomers), are 
useful features to determine the different behavior of the heat capacity as a function of the cluster 
size. The effect of the range of the interatomic forces is studied by comparing the heat capacities of 
the Nass and Cdss clusters. It is shown that by decreasing the range of the many-body interaction, 
the distribution of isomers characterizing the PEL is modified appropriately to generate a negative 
heat capacity in the Cdss cluster. 
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I. INTRODUCTION 

Negative microcanonical heat capacity in atomic and 
molecular clusters was theoretically predicted by consid- 
ering simple models of the distribution of local minima 
that characterize the potential energy landscape (PEL) 
of clusters^. In that study, it was found that for high 
values of the parameter involving the ratios of the vibra- 
tional frequencies corresponding to the global and local 
isomers, the caloric curve displays an S-shaped loop, with 
a negative heat capacity in the vicinity of the melting 
point 1 . In another study on the solid-liquid transition 
of clusters^, it was shown that in microcanonical sim- 
ulations of Lennard- Jones clusters, an increase in total 
energy causes a temperature reduction. This effect was 
related to the broadening of the cluster kinetic-energy 
distribution toward lower energy values^. 

Although the existence of negative heat capacity in 
physical systems like stars or star clusters^, and in 
fragmenting nuclei^ is well-known, this peculiar effect 
gained a lot of interest in the field of atomic and molec- 
ular clusters due to recent experimental results where 
a negative heat capacity was measured for a 147-atom 
sodium cluster—. In this study, the photofragmentation 
mass spectra was used to measure the internal energy 
of free, mass selected clusters with known temperature. 
These measurements were used to determine the micro- 
canonical caloric curve of the Na + i47 that shows the char- 
acteristic S-shaped (backbending) feature, indicating a 
negative heat capacity^. The negative value of the mi- 
crocanonical heat capacity was interpreted by consider- 
ing that a finite system upon melting tries to avoid partly 
molten states and prefers to convert some of its kinetic 
energy into potential energy^. This peculiarity has been 
attributed to the non-additivity of the total energy of a 
cluster with finite sizeZA 

Microcanonical heat capacities of metal clusters 
have been theoretically investigated using constant- 



energy molecular dynamics (MD) with many-body 
potentialsSiifliii and an orbital-free version of the first- 
principles MD methodiSii^. In these studies, heat capac- 
ities of fee transition and noble metal clusters with up to 
23 atoms were calculated to characterize their melting- 
like transition 9 . In another study, on the melting of 
sodium clusters^, the microcanonical caloric curve of 
the Na55 cluster was obtained. However, in not one of 
these studies was a signature of a negative heat capac- 
ity found. Similar results, indicating the non-existence of 
a negative heat capacity in constant-energy orbital-free 
first-principles MD simulations of larger sodium clusters 
(Na55, Nag2, and Nai42), were obtained^. Nevertheless, 
in such calculations the simulation time employed was 
too short to obtain converged results— On the other 
hand, in microcanonical MD simulations of Aljy, N = 7, 
13, 55, and 147, clusters, a negative heat capacity was 
obtained for the larger AI147 cluster—. 

In the present work, motivated by the availability of ex- 
perimental techniques that allow the measurement of the 
microcanonical heat capacity and other thermal proper- 
ties of mass selected metal clustersi^i^i^, we theoret- 
ically investigate the behavior of the heat capacity of 
sodium clusters in the size range of 13-147 atoms. In 
our approach, constant-energy MD simulations are per- 
formed using a phenomenological many-body potential 
that mimics the metallic bonding of sodium clusters. 
This approximation allows us to use simulation times 
of the order of ~ 50 ns, in order to obtain converged 
averages of the microcanonical heat capacity and other 
cluster thermal properties. Our main objective is to gain 
additional insights into the conditions that determine if 
a cluster has a negative heat capacity. The main find- 
ing of this work shows that the width of the distribution 
function of the kinetic energy and the spread of the dis- 
tribution of the potential energy minima (isomers), char- 
acterizing the PEL, are useful features to determine the 
signature of the cluster heat capacity. In section II, wc 
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provide the theoretical background on which this study 
is based. The results and their discussion are given in 
section III, and section IV contains a summary and the 
conclusions of this work. 



II. THEORETICAL BACKGROUND 

The heat capacity and temperature of sodium clusters 
as a function of the cluster total energy are calculated 
through constant-energy MD simulations u sing the mi- 
crocanonical expressions derived in Refs. 1171113 : 
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,4=0.01595 eV, f=0. 29113 eV, r =6.99 bohr, p=10.13, 
and q=l. 3Q2£. This phenomelogical many-body poten- 
tial has been used to study the melting-like transition 
in sodium clusters of different sizes using Monte Carlo 
(MC)iS and constant-energy MD simulations^. A good 
qualitative agreement has been obtained between struc- 
tural and thermal properties calculated using the Gupta 
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and those generated from first-principles 



methodsiSii 3 .. An additional advantage in using this po- 
tential is that it allows simulation times of the order of 
50 ns, necessary to obtain fully converged time averages 
in the melting-like transition region. 



III. RESULTS AND DISCUSSION 
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where K is the kinetic energy of the cluster, fee is the 
Boltzmann constant, and (...) denotes a time average. 
In these formulas, 3iV was changed to 3A-6, the num- 
ber of degrees of freedom of the system, since the cal- 
culations are performed for a non-translating and non- 
rotating cluster in a three-dimensional space (the posi- 
tion of the center of mass was fixed and the total mo- 
mentum was held to zero during the simulations). 

In our implementation of the constant-energy MD 
method, the Newton's equations of motion are solved 
with the Verlet algorithmic using a time step of 2.4 fs, 
which provides total energy conservation within 0.001 
%. A typical calculation consists in heating up a clus- 
ter from its lowest-energy solid-like configuration until 
it transforms into a liquid-like cluster. To simulate this 
procedure the cluster total energy is increased in a step- 
like manner by scaling up the atomic velocities. For each 
initial condition the cluster was equilibrated during 10 4 
time steps and the time averages of the physical quanti- 
ties are calculated using 10 7 time steps. This averaging 
time is increased by a factor of 2 when the cluster is in the 
region of the solid-to-liquid transition in order to ensure 
the calculation of fully converged averages. 

To model the metallic bonding in sodium clusters we 
used the many-body Gupta potential^, which is based 
on the second moment approximation of a tight-binding 
hamiltonian^i. Its analytical expression is given by: 
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where r , A, £, p, and q are adjustable parameters 2 ^. 
For sodium clusters these parameters have been fit- 
ted to band structure calculations 22 . Their values are: 



The microcanonical heat capacities of the Na^, N = 
13, 20, 55, 135, 142, and 147, clusters, calculated using 
Eq. (1), are displayed in Fig. □ For N = 13, 20, and 
55, they are continuous functions of the cluster total en- 
ergy showing a maximum value that is characteristic of 
a melting-like transitioni2ii2iiii2ii. On the other hand, 
the heat capacity of the larger clusters (N = 135, 142, 
and 147) shows two discontinuity points and a continuous 
negative-valued interval between them. This peculiar be- 
havior in the heat capacity is related with a backbending 
loop in the caloric curve (temperature as a function of the 
total energy )ii2i2i 2 ^. In fact, our calculated microcanon- 
ical caloric curves of Nai 35 , Nai 42 , and Nai.47 show the 
backbending loop at the same energies where the heat ca- 
pacity takes negative values (see the caloric curves shown 
in Fig. 2 of Ref. l23j) . In previous studies, the negative 
slope (backbending loop) of the microcanonical caloric 
curve has been attributed to a peculiar behavior of the 
cluster entropy as a function of energy that shows a dent 
with inverted curvature in the region of the solid-liquid 
transition 1 ! 2 ! 7 ! 8 ! 25 . 

In the present work, we analyze the behavior of the 
microcanonical heat capacity of sodium clusters from a 
different perspective. First, we consider Eq. (1) and 
obtain the condition to have a negative value in the heat 
capacity: 
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Figure [21 display the values of Ze (black dots) as a func- 
tion of the cluster total energy E that were calculated 
from a time average, using the MD trajectories. In the 
same scale the threshold value Z c = (3N — 6)/(3A — 8) 
for each cluster size is given. In Fig. [5], it can be graph- 
ically seen how the relative difference between Ze and 
Z c changes with the cluster size. For the three smaller 
clusters (see panels (a), (b), and (c) in Fig. [21 ) the Z E 
values do not overcome the threshold value, whereas for 
the three larger clusters there is a range of total energy 
where Ze satisfy the condition to have negative heat ca- 
pacities (see panels (d), (e), and (f) in Fig. 
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In order to investigate what determines a negative 
value of the heat capacity, we consider the quantity Ze- 
This is the product of the averages of the kinetic energy 
and of the inverse of this quantity, and therefore, its value 
will depend on the distribution function of the kinetic en- 
ergy, gE(K). The average of any function of the kinetic 
energy f(K) can be obtained through the following ex- 
pression: 

(f E ) g = J f{K)g E {K)dK. (6) 

Since the distribution function of the kinetic energy, 
gE{K), determines the behavior of Ze, it is useful to 
analyze gE{K) at different values of the cluster total en- 
ergy E. The calculation of this quantity is straightfor- 
ward from the constant-energy MD simulations. Figure 
13 shows as a function of the normalized mean devia- 
tion 8 K — (K - (K))/ (K) for three different energies, 
corresponding to the cases where the cluster is in the 
solid- and liquid-like phases, and at the middle of the 
melting-like transition. The analysis of gs as a function 
of SK, instead of K, has the advantage that it allows 
the comparison, on the same scale, of the lineshapes of 
this function at different cluster energies and for different 
cluster sizes. 

As a general trend, it is found that gE becomes nar- 
rower for increasing cluster sizes, indicating a larger rela- 
tive dispersion of the kinetic energy values for the smaller 
clusters. This result is expected since it confirms the in- 
crement of fluctuations in kinetic energy of a physical 
system that decrease in size. A common characteristic of 
gE, existing in the six clusters investigated, is the larger 
broadening of the distribution function when the clus- 
ter is at the melting-like transition. At lower (solid-like 
phase) and higher (liquid- like phase) energies, the width 
of gE is smaller, whereas at the phase transition the fluc- 
tuations in kinetic energy, as expected, should increase. 

For the three smaller clusters which do not have neg- 
ative heat capacity, gE^SK) shows a nearly symmetric 
lineshape independent of the cluster energy (see panels 
(a), (b), and (c) of Fig. |3|). In contrast, Nai42 and Nai47, 
that show a negative heat capacity, have a distribution 
function gE{SK) with a shoulder towards positive values 
of SK, at energies in the middle of the melting region (see 
panels (e) and (f) of Fig. [21). Although this difference in 
the distribution function of the kinetic energy could be 
a useful feature to determine the existence of a negative 
heat capacity, the Nai35 cluster would be an exception 
to this rule since its gs(SK) does not show a resolved 
shoulder in its lineshape (see panel (d) in Fig. 0), but it 
has negative heat capacity. 

On the other hand, a characteristic of gE that would 
be useful to determine the sign of the heat capacity is the 
width of the distribution function, which can be obtained 
through its second moment: 

({8Kf) g = J (5K) 2 g E (5K)d(5K). (7) 



The second moment of gE corresponds to the second term 
in the expansion of Ze, which can be obtained from the 
left hand term of Eq. (5)iZ»i2c 

Z E = l + ((SK) 2 ) g + .... (8) 

By taking terms up to second order in this expansion, 
assuming that (SK) « 1, Ze can be approximated by: 

Z E ^zf = \ + {{5Kf) g . (9) 

(2) 

Since the Z E values can be calculated using gE and Eq. 
(7), it is possible to check the validity of this approxima- 
tion, which can only be applied to systems with a finite 
number of particles. Fig. [3 shows the values of Z E ^ 
(stars) as a function of the cluster energy. It can be seen 
that the difference with Ze is small for the three smaller 
clusters and negligible for the three larger ones. Then, 

(2) 

Z E can be considered as a quantitative measure of the 
width (second moment) of the distribution function of the 
kinetic energy, and can be used to determine the sign of 
the heat capacity. Figure ^ shows the maximum values 

(2) 

of Z E (black dots) calculated for energies at the mid- 
dle of the melting transition, and their comparison with 
the threshold values Z c as a function of the cluster size 
(full line). From this figure, it is obvious that although 
the width of gE is relatively large for the three smaller 
clusters, the corresponding maximum values of Z E are 
below the threshold Z c , and therefore these clusters do 
not show a negative heat capacity. On the other hand, 
the width of gE for Nai35, Nai42, and Nai47 is smaller 
than for Nai3, Na2o, and Na55, however, Z c is a faster 
decreasing function of the cluster size N , such that the 

(2) . 

maximum of Z E lies above the threshold values, indi- 
cating that the larger clusters have negative heat capac- 
ities. Then, the above results suggest that the width of 
the distribution function of the kinetic energy is a useful 
property to determine the sign of the heat capacity of 
clusters. However, this quantity not only depends on the 
cluster size, but also on the characteristics of the PEL. 

To illustrate the importance of the topology of the PEL 
we have investigated the behavior of the heat capacity of 
55-atom clusters using the many-body Gupta potential, 
shown in Eqs. (3) and (4), for the different metals listed 
in Tables I and III of Ref. 21, and in Table II of Ref. 
I22I Our results show that the Cdss cluster (with the fol- 
lowing parameter values^: ,4=0.0416 eV, £=0.4720 eV, 
p=13.639, and q=3.908) has a negative heat capacity, 
but none of the other 55-atom clusters show this pecu- 
liarity. The upper insets of Fig. 0] show the calculated 
caloric curve with the corresponding backbending loop 
and the heat capacity with negative values for a range of 
total energy values of the Cdss cluster. The inset at the 
lower right corner of Fig. 4 shows gE as a function of 
the normalized mean deviation of the kinetic energy for 
both, the Na55 and Cdss, clusters. It can be seen that 
the broadening of gE is larger in Cdss than in Nass, such 
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that, it generates a maximum value of Z E (this value 
corresponds to the point represented by a star in Fig. 4) 
that overcome the threshold Z c , and consequently, the 
Cdss cluster display a negative heat capacity. This com- 
parison with the Nass cluster which does not show this 
peculiarity, indicates that although both clusters have 
the same size, their dynamical properties defined by their 
corresponding PEL's, generate different behavior in their 
heat capacities. 

In order to investigate the influence of the PEL on the 
different widths of the distribution function of the kinetic 
energy of the Nass and Cdss clusters, further studies are 
necessary. In this direction, the calculation of short-time 
averages of the kinetic energy and periodic quenchings of 
instantaneous configurations during the MD trajectories 
allow us to obtain the distribution of potential energy 
minima (isomers) that are accessible at different clus- 
ter energies^. Figure [S] display the normalized distri- 
bution of potential energy minima, obtained by period- 
ical quenchings using MD trajectories at a total energy 
where the cluster is at the middle of the melting transi- 
tion, for the Nass and Cdss clusters. It can be notice that 
the number of isomers with higher energy relative to the 
global minimum, are larger for the Cdss cluster in com- 
parison with the results obtained for Nass. This result 
can be explained by taking into account that the range 
of the interatomic forces is shorter in Cd than in Na clus- 
ters, mainly due to the higher value of the q parameter in 
the many-body Gupta potential^. The physical reason 
for the larger number of minima at short range is the loss 
of accessible configuration space as the potential wells be- 
come narrower, thus producing barriers where there arc 
none at long ranged. 

To show how the distribution of potential energy min- 
ima determine the broadening of the distribution func- 
tion of the kinetic energy, we approximate the complex 
topology of the PEL by a set of independent harmonic 
potential wells in the 3N — 6 dimensional space. Each 
one of these wells is associated to the different potential 
energy minima forming the distribution of isomers shown 
in Fig. 5. For each potential energy minimum denoted 
by I, the distribution function of the kinetic energy at a 
total energy E, in the harmonic approximation, is given 
by^±£: 

g E ,i(K) = Ci(E — Ai — K) 1 ^ 1 K^ 1 : (10) 

where A; is the potential energy of the isomer I, relative 
to the potential energy value of the lowest-energy isomer, 
and Ci is a normalization constant such that: 

/ g E .i(K)dK = 1. (11) 
Jo 

The distribution function of the kinetic energy gE.har, at 
a total energy E, corresponding to the whole PEL can be 
constructed by adding up the contribution of each har- 
monic potential well, weighted by the probability, ug j, 
of finding a given isomer during the quenching from the 



MD trajectories. This probability is given by the height 
of the distribution shown in Fig. 5. Then, gg h ar is given 
by: 

9E,har = X! U E,l9E,l(K), (12) 
1=1 

with 

£"wb,1 = 1. (13) 
1=1 

By using the data from the whole distribution of isomers 
in Fig. 5, the distribution function of the kinetic energy 
9E,har was calculated for the Cdss and Nass clusters. 
They are displayed in the insets of Fig. 5 (full lines). A 
comparison between 9E,har 

and the exact gE (obtained 
from the MD simulation and displayed in the lower right 
inset of Fig. 4) shows that there is a good agreement be- 
tween the two distribution functions. This indicates that 
gE is determined mainly from the number of isomers and 
the probability to found them (height of the distribu- 
tion), rather than from other features of the PEL like 
saddle points. The advantage in introducing gE,har m 
this discussion is related with the fact that it is possible 
to analyze the broadening of this distribution function 
of the kinetic energy by considering different subsets of 
potential energy minima (isomers). This is useful to de- 
termine what regions of the PEL are more relevant to 
increase the width of gE,har, and investigate the appear- 
ance of the negative heat capacity. The insets of Fig. 5 
show three partial distribution functions gE,har, consid- 
ering different subsets of isomers corresponding to three 
intervals of low (L), medium (M) and high (H) potential 
energy values. By analyzing the relative contribution of 
these subsets to the width of gE,har for the Cdss and 
Nass clusters, it is found that the larger broadening in 
the cadmium cluster is mainly due to the larger contribu- 
tion of the isomers in the range of high potential energy 
which are spreaded along a larger interval than those cor- 
responding to the Nass cluster. The width of gE,har for 
the Nass cluster is smaller since there are proportion- 
ally less isomers with high potential energy, and they are 
extended over a shorter interval of values. As was men- 
tioned above, the physical reason for this difference in the 
distribution of isomers between the Cdss and Nass clus- 
ters is the shorter range of the many-body forces existing 
in the cadmium cluster as compared with those present 
in the sodium cluster. A similar result was obtained for 
55-atom clusters using a pairwise Morse potential for dif- 
ferent values of the range of the interatomic forces^. In 
that case the backbending loop in the caloric curve (neg- 
ative heat capacity) was obtained using a Morse poten- 
tial with a range of the interatomic forces that is shorter 
than the range characteristic of alkali metals which have 
long-ranged interactions^. 

Therefore, if a detailed characterization of the distribu- 
tion of isomers forming the PEL of clusters is performed, 
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the broadening of qe may be estimated, and by the com- 

(2) 

parison of the corresponding Z E and Z c values, it would 
be possible to predict the sign of the heat capacity of 
clusters. 



IV. SUMMARY 

The microcanonical heat capacity of sodium clusters 
has been calculated using constant-energy MD simula- 
tions and the many-body Gupta potential. Negative val- 
ues for the heat capacity at energies where the cluster 
is at the melting-like transition were found for Nai35, 
Nai42, and Nai47. The smaller sodium clusters Na^r, 
iV=13, 20, and 55, do not show this peculiarity. An 
analysis of the calculated distribution function of the ki- 
netic energy gE for the six clusters investigated, shows 
that the width of this distribution function is a useful 
feature to determine the sign of the heat capacity. It 
was found, that although the broadening of gE is larger 
for the smaller clusters, it is not enough to overcome the 
corresponding threshold value to obtain a negative heat 
capacity. However, since this threshold is a fast decreas- 



ing function of the cluster size, the broadening of gE in 
the larger clusters is enough to generate a negative heat 
capacity. 

It was also shown that the broadening of gE depends 
on the distribution of potential energy minima that char- 
acterize the PEL of clusters. Specifically, as the range of 
the many-body interactions is decreased (like in the case 
of Cd clusters), the number of local minima with higher 
energy increases generating a larger broadening in gE, 
and consequently a negative heat capacity. The analysis 
presented in this paper shows how the complex topology 
of the PEL can be explored to extract the main features 
that determine the sign of the heat capacity of metal 
clusters. 
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FIG. 1: Heat capacity of Na N , N = 13 (a); 20 (b); 55 (c); 135 
(d); 142 (e); and 147 (f) clusters. The cluster energy is cal- 
culated taking as reference the value of the binding energy of 
the most-stable (lowest-energy) configuration given in Table 
I of Ref. 



FIG. 3: Distribution function of the kinetic energy for Najv, 
N = 13 (a); 20 (b); 55 (c); 135 (d); 142 (e); and 147 (f) 
clusters. The three curves displayed in each panel correspond 
to the solid- (lower energy), melting- (intermediate energy), 
and liquid-like (higher energy) phases. 



FIG. 2: Energy dependence of the Ze (black dots) and Z E 
(stars) values for Najv, N = 13 (a); 20 (b); 55 (c); 135 (d); 142 
(e); and 147 (f) clusters. The Ze values were calculated using 
Eq. (5) whereas the Z% values, which are an approximation 
of Ze according to Eq. (9), were obtained using the second 
moment of the distribution function of the kinetic energy. See 
the related text for an explanation of the difference between 
these quantities. The cluster energy is calculated taking as 
reference the value of the binding energy of the most-stable 
(lowest-energy) configuration given in Table I of Ref. l23i 
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FIG. 4: Comparison of the maximum values of Z E (black 
dots for the Naw clusters) and the threshold Z c (continuous 
line) as a function of the cluster size. The star shows the 
maximum value of Z E ' for the Cdss cluster. The upper in- 
sets show the energy dependence of the caloric curve and the 
heat capacity of the Cdss cluster. The inset at the low right 
corner shows the §e as a function of the normalized mean 
deviation of the kinetic energy of the Cdss cluster, and its 
comparison with Nass, calculated at E = 4.37 eV and E = 
2.71 eV, respectively. 



FIG. 5: Normalized distribution of potential energy min- 
ima for the Nass and Cdss clusters. This distribution was 
obtained from two thousand quenchings separated by 10000 
time steps during a MD trajectory of 20 million time steps. 
The values of E=4.37 eV (Cd 55 ) and E=2.71 eV (Na 55 ) cor- 
respond to the cluster energies when they are at the middle 
of the melting-like transition. The vertical dashed lines sepa- 
rate intervals of low (L), medium (M), and high (H) potential 
energy. The insets show the total (full line), including the L, 
M, and H intervals, and partial (taking different subsets of 
isomers) gE,har, as a function of the normalized mean kinetic 
energy. 
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